;-------------------------------------------------------------LICENSE--------------------------------------------------------------;
;                                                                                                                                  ;
;The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    ;
;and Message Passing Interface (MPI) parallelization.                                                                              ;
;                                                                                                                                  ;
;Copyright (C) 2012                                                                                                                ;
;Ronglin Jiang                                                                                                                     ;
;rljiang@ssc.net.cn                                                                                                                ;
;585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             ;
;                                                                                                                                  ;
;This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         ;
;as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             ;
;                                                                                                                                  ;
;This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    ;
;MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        ;
;                                                                                                                                  ;
;You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     ;
;Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   ;
;                                                                                                                                  ;
;-------------------------------------------------------------LICENSE--------------------------------------------------------------;

function current_density,bx,by,bz,dx,dy,dz,nx,ny,nz

on_error,2

current=dblarr(nx,ny,nz,3)

current[*,*,*,0]=((shift(bz,0,-1,0)-shift(bz,0,1,0))/repli3d(dy,dir='y',nx=nx,ny=ny,nz=nz)-(shift(by,0,0,-1)-shift(by,0,0,1))/repli3d(dz,dir='z',nx=nx,ny=ny,nz=nz))/4/!pi
current[*,*,*,1]=((shift(bx,0,0,-1)-shift(bx,0,0,1))/repli3d(dz,dir='z',nx=nx,ny=ny,nz=nz)-(shift(bz,-1,0,0)-shift(bz,1,0,0))/repli3d(dx,dir='x',nx=nx,ny=ny,nz=nz))/4/!pi
current[*,*,*,2]=((shift(by,-1,0,0)-shift(by,1,0,0))/repli3d(dx,dir='x',nx=nx,ny=ny,nz=nz)-(shift(bx,0,-1,0)-shift(bx,0,1,0))/repli3d(dy,dir='y',nx=nx,ny=ny,nz=nz))/4/!pi

current[0,1:ny-2,1:nz-2]=current[1,1:ny-2,1:nz-2]
current[nx-1,1:ny-2,1:nz-2]=current[nx-2,1:ny-2,1:nz-2]
current[*,0,1:nz-2]=current[*,1,1:nz-2]
current[*,ny-1,1:nz-2]=current[*,ny-2,1:nz-2]
current[*,*,0]=current[*,*,1]
current[*,*,nz-1]=current[*,*,nz-2]

return,current

end